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Limit possibilities of observer-based synclironization systems under information constraints (lim- 
ited information capacity of the coupling channel) are evaluated. We give theoretical analysis for 
multi-dimensional drive-response systems represented in the Lurie form (linear part plus nonlin- 
earity depending only on measurable outputs). It is shown that the upper bound of the limit 
synchronization error (LSE) is proportional to the upper bound of the transmission error. As a 
consequence, the upper and lower bounds of LSE are proportional to the maximum rate of the 
coupling signal and inversely proportional to the information transmission rate (channel capacity). 
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applied to synchronization of two chaotic Chua systems coupled via a channel with limited capacity. 
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I. INTRODUCTION 



Chaotic synchronization has attracted the attention of 
researchers since the 1980s and is still an area of 

active research IS ■ Recently information-theoretic 
concepts were ap plie d to analyze and quantify synchro- 
nization 0, H IE IS HI' In H, Hi mutual information 
measures were introduced for evaluating the degree of 
chaotic synchronization. In 0, 0| the methods of sym- 
bolic dynamics were used to relate synchronization pre- 
cision to capacity of the information channel and to the 
entropy of the drive system. Baptista and Kurths ^Jl 
introduced the concept of a chaotic channel as a medium 
formed by a network of chaotic systems that enables in- 
formation from a source to pass from one system (trans- 
mitter) to another system (receiver). They character- 
ized a chaotic channel by the mutual information (dif- 
ference between the sum of the positive Lyapunov ex- 
ponents corresponding to the synchronization manifold 
and the sum of positive exponents corresponding to the 
transverse manifold). However, in existing papers limit 
possibilities for the precision of controlled synchroniza- 
tion have not been analyzed. 

Recently the limitations of control under constraints 
imposed by a finite capacity information channel have 
been investigated in detail in the control theoretic litera- 
ture, see ^l,|l3l[i3jliEI3 and references therein. It was 
shown that stabilization under information constraints is 
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possible if and only if the capacity of the information 
channel exceeds the entropy production of the system at 
the equilibrium [Hill Eg. In [Hill a general state- 
ment was proposed, claiming that the difference between 
the entropies of the open loop and the closed loop sys- 
tems cannot exceed the information introduced by the 
controller, including the transmission rate of the infor- 
mation channel. However, results of the previous works 
on control system analysis under information constraints 
do not apply to synchronization systems since in a syn- 
chronization problem trajectories in the phase space con- 
verge to a set (a manifold) rather than to a point, i.e. the 
problem cannot be reduced to simple stabilization. 

In this paper we establish limit possibilities of observer- 
based synchronization systems under information con- 
straints. Observer-based synchronization systems are 
used when only one phase variable is available for mea- 
surement and coupling. Such systems are well studied 
without information constraints [l^ l20l |2]| . Here we 
present a theoretical analysis for n-dimensional drive- 
response systems represented in the so called Lurie form 
(linear part plus nonlinearity, depending only on measur- 
able outputs) . It is shown that the upper bound of the 
limit synchronization error (LSE) is proportional to the 
upper bound of the transmission error. As a consequence, 
the upper and lower bounds of LSE are proportional to 
the maximum rate of the coupling signal and inversely 
proportional to the information transmission rate (chan- 
nel capacity). Optimality of the binary coding for coders 
with one-step memory is established. 

Note also that it was claimed in some papers that if the 
capacity of the channel is larger than the Kolmogorov- 
Sinai entropy of the driving system, then the synchro- 
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nization error can be made arbitrarily small. Such a 
claim is based upon the noisy channel theorem of Shan- 
non information theory stating that, if the source entropy 
is smaller than the channel capacity, then the data gen- 
erated by the source can be transmitted over the chan- 
nel with negligible probability of error. However, it is 
known that to transmit data with a sufhciently small 
error a sufficiently long codeword and a long transmis- 
sion time is needed. During such a long time an unstable 
chaotic trajectory may go far from the its predicted value 
and synchronization may fail. Therefore analysis of the 
system precision under information constraints requires 
more subtle arguments which are provided in this paper 
based on Lyapunov functions and coding analysis. 



II. DESCRIPTION OF OBSERVED-BASED 
SYNCHRONIZATION SYSTEM 

To simplify exposition we will consider unidirectionally 
coupled systems in the so-called Luric form: right-hand 
sides are split into a linear part and a nonlincarity vector 
depending only on the measured output. Then the drive 
system is modelled as follows: 



X = Ax + Lp{y), y = Cx, 



(1) 



where x is an n-dimensional (column) vector of state vari- 
ables, y is the scalar output (coupling) variable, A is an 
n X n-matrix, C is n x 1 (row) matrix, (p{y) is a contin- 
uous nonlincarity. We assume that the vector of initial 
conditions xq = x{0) belongs to a bounded set such 
that all the trajectories of the system ^ starting in 
are bounded. Such an assumption is typical for chaotic 
systems. 

The response system is described as a nonlinear ob- 
server 



X = Ax + (p{y) + L{y - y), y = Cx, 



(2) 



where L is the vector of the observer parameters (gain). 
Apparently, the dynamics of the state error vector e(t) = 
x{t) — x{t) is described by a linear equation 



A 



Le, y 



Cx, 



(3) 



where Al = A - LC . 

As is known from control theory, e.g. |22l |. 
if the pair {A, C) is observable, i.e. if 
rank(C^,A^C^,...,(A^)"-iC^) = n, then there 
exists L providing the matrix A^ with any given 
eigenvalues. Particularly, all eigenvalues of Al can have 
negative real parts, i.e. the system Q can be made 
asymptotically stable and e{t) Q ast ^ oo. Therefore, 
in the absence of measurement and transmission errors 
the synchronization error decays to zero. 

Now let us take into account transmission errors. We 
assume that the observation signal y{t) is coded with 
symbols from a finite alphabet at discrete sampling time 



instants tk = kTs,k = 0,1,2,..., where is the sam- 
pling time. Let the coded symbol yk = y{tk) be trans- 
mitted over a digital communication channel with a fi- 
nite capacity. To simplify the analysis, we assume that 
the observations are not corrupted by observation noise; 
transmissions delay and transmission channel distortions 
may be neglected and the coded symbols are available at 
the receiver side at the same sampling instant tk = kTs- 
Assume that zero- order extrapolation is used to convert 
the digital sequence yk to the continuous-time input of 
the response system y{t), namely, that yit) — yk as kTs < 
t < (k + l)Ts. Then the transmission error is defined as 
follows: 



Syit) = y{t)-m- 



(4) 



In presence of the transmission error, equation (0) 
takes the form 



^Ae + (^(y) ~ ^{y + Sy{t)) - LySy{t) 



(5) 



Our goal is to evaluate limitations imposed on the syn- 
chronization precision by limited transmission rate. To 
this end introduce an upper bound of the limit synchro- 
nization error Q — sup lim ||e(i)||, where e{t) is from jSJl, 

t — ^oo 

II • II denotes the Euclidian norm of a vector, and the 
supremum is taken over all admissible transmission er- 
rors. In the next two sections we describe encoding and 
decoding procedures and evaluate the set of admissible 
transmission errors Sy{t) for the optimal choice of coder 
parameters. It will be shown that 6y{t) is bounded and 
does not tend to zero. 



III. CODING PROCEDURES 

At first, consider the memoryless (static) encoder with 
uniform discretization and constant range. For a given 
real number M > and positive integer v € Z define a 
uniform scaled coder to be a discretized map q^^M '■ R — *■ 
M as follows. Introduce the range interval T = [—M,AI] 
of length 2M and the discretization interval of length 
S = 2^~''M and define the coder function qu,M[y) as 



<1u,m{v) 



6-{S-^y), if \y\<M, 
M sign(y), otherwise, 



(6) 



where (•) denotes round-up to the nearest integer func- 
tion, sign(-) is the signum function: sign(y) = 1, if y > 0, 
sign(y) = -1, if y < 0. Evidently, \y-qi,,M{y)\ < S/2 for 
all y such that y : |y| < M + 6/2 and all values of qv,M{y) 
belong to the range interval X. Notice that the interval I 
is equally split into 2" parts. Therefore, the cardinality 
of the mapping q^^M image is equal to 2^-1-1, and each 
codeword symbol contains R — log2(2'' -f 1) bits of infor- 
mation. Thus, the discretized output of the considered 
encoder is found as y = qv.kiiy)- We assume that the 
encoder and decoder make decisions based on the same 
information. 
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In a number of papers more sophisticated enco ding 
schemes have been proposed and analyzed, see 0, 

I25I I for example. The underlying idea for coders of 
this kind is to reduce the range parameter M, replac- 
ing the symmetric range interval I by the interval yt+i, 
covering some area around the predicted value for the 
{k + l)th observation yk+i, Vk+i G yk+i- If the length 
of yk+i is small compared with the full range of possible 
measured output values y, then there is an opportunity 
to reduce the range parameter M and, consequently, to 
decrease the coding interval S preserving the bit-rate of 
transmission. To realize this scheme, memory should be 
introduced into the encoder. Using such a "zooming" 
strategy it is possible to increase coder accuracy in the 
steady-state mode, and, at the same time, to prevent 
coder saturation at the beginning of the process. 

In this paper we use a simple version of such an encoder 
having one-step memory and time-based zooming. To 
describe it we introduce the sequence of central numbers 
Ck, k — 0,1,2, . . . with initial condition cq — 0. At step 
k the encoder compares the current measured output yk 
with the number Ck, forming the deviation signal dyk = 
yk ~ Ck- Then this signal is discretized with a given v 
and M = Mk according to ©. The output signal 



dyk = qy,Mk{dyk) 



(7) 



is represented as an i?-bit information symbol from the 
coding alphabet and transmitted over the communication 
channel to the decoder. Then the central number Ck+i 
and the range parameter Mk are renewed based on the 
available information about the driving system dynamics. 
We use the following update algorithms: 



Ck+i^Ck+dyk, Co = 0, A: = 0,1, 



Mfe^CMo-MooV' + Moo, fc = 0,l,. 



(8) 



(9) 



where < p < 1 is the decay parameter, M^o stands for 
the limit value of Mk- The initial value Mq should be 
large enough to capture all the region of possible initial 
values of y^. 

The equations ®, iQ, describe the encoder al- 
gorithm. The same algorithm is realized by the de- 
coder. Namely, the decoder calculates the variables Ck, 
Mk based on received codeword flow similarly to Ck, Mk- 



IV. CODER OPTIMIZATION 

We now find a relation between the transmission rate 
and the achievable accuracy of the coder-decoder pair, 
assuming that the growth rate of y{t) is uniformly 
bounded. Obviously, the exact bound Ly for the rate 
of y{t) is Ly — supjCil, where x is from To an- 

alyze the coder-encoder accuracy, evaluate the upper 

bound A = sup of the transmission error 5y{t) — 

t 



y{t) — y(t). Consider the sampling interval [tk,tk+i]- It 

is clear that \5y{tk)\ does not exceed 5/2. Additionally, 

the error may increase from tk to tk+i due to change 

of y{t) by a value not exceeding sup \y{t) — y{tk)\ 

tk <t<tfc+i 

< / \y{'r)\dT < J Lydr = LyTs. Therefore the total 

tk tk 

transmission error for each interval [tk, ifc+i] satisfies the 
inequality: 



\Sy{t)\ < d/2 + LyTs 



(10) 



Inequality (|10|l shows that in order to meet the inequality 
|(5y(t)| < A for all t, the sampling interval Tg should 
satisfy condition 



T, < A/Ly. 



(11) 



Furthermore, if the condition Hll|) holds, the given 
bound for the coding error will be guaranteed if the 
coding interval S is appropriately chosen, namely, 5 < 
2 A — 2LyTs. It provides the lower bound for the trans- 
mission bit-per-step rate B = log2(2Af/(5 +1): its value 

should not be less than log2 ( a + 1 ) . Therefore, 



the coder with range 2M, coding interval 6 and sampling 
period Tg ensures the total transmission error A if Hll() 
holds and the transmission rate satisfies inequality 



R > log2 



M 



A - LyTs 



(12) 



It follows from (|11() . H12() that if Tg is sufficiently small 
and R is sufficiently large, then an arbitrarily small value 
of A can be assured. 

Let us now optimize the coder parameters to achieve 
the minimum bound for the error A. As seen from H12|l . 
reduction in the coder range 2M results in a reduction 
in the transmission rate R and channel capacity R* . On 
the other hand, to prevent coder saturation, M should 

not be less than sup\Sy{tk)\ — S/2 = A — S/2. Taking 

kei, 

into account that S — 2^~'^M, we arrive at the following 
formula for the minimal admissible range: 



M = 



2^ + 1 



-A. 



(13) 



Remark 1. At the initial stage of the system evolution 
the error \5y \ may exceed the bound A, because the initial 
value y(0) is not known. This leads to the transient mode 
of the system behavior. The zooming strategy may be 
efficient at this stage, providing the following recipe for 
the choice of the coder parameters: Mq — My, Mryo — 
2''A/{2'' + 1), where My = sup \y{t)\. 

Now optimize the coder w.r.t. T^. Consider the steady- 
state mode when \y{t) — Ck\ < A for each time interval 
i € [tk:tk+i)- Let M be found from 113|) . Introduce a 
real number e as e = LyT^j A, (evidently, < e < 1) and 
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rewrite the lower bound R* for R in the form 

R* = log ( %- + 1 ) . (14) 

Defining the bit-per- second rate R = R/T^ and its lower 
bound R* we have from H14|l : 

= ^ (2.^ + 1)(1 -.) +')' 

Now the optimization of the coder is reduced to 
the following minimization problem: Find {e*,i^*) = 
arg min R{£,v). Since the right-hand side of (I15II is 

ee(0,l) 

strictly growing in the optimal value of is i^* = 0. 
This means that the binary coding scheme gives the op- 
timal transmission rate R* — 1 bit per step, which yields 
M* = A/2 as an optimal value for M, and the signum 
function as an optimal coder function: y = y ^^S^V- 

For the optimal value of v we have R* — {Ly/ A)r{e), 
where 

K.)4log.(^ + l). (16) 
Let r* — min r(e). It is easy to see that this mim- 

0<e<l ^ ^ 

imum exists and satisfies the transcendental equation 
dr{e)/de = 0. Numerical one-dimensional minimization 
yields r* = r{e*) « 1.688, where e* « 0.5923. 
Therefore, the optimal sampling time T* is 

T;=e*A. (17) 

Then the minimal channel bit-rate R* = l/Tg is 

R*^r*^, (18) 

and this bound is tight for the considered class of coders. 
The relation H18|l can be rewritten as 

RA > r*Ly, (19) 

playing the role of an uncertainty relation between the 
propagation rate of information and the transmission er- 
ror. 

Remark 2. The limits of synchronization error may 
be different if a more sophisticated coder is used, e.g. a 
first-order coder with linear extrapolation of the signal 
or nth order coder with predictive model of the drive 
system. For example, if a full order observer is admit- 
ted at the transmitter side and there are n channels 
for simultaneous transmission of the n-dimensional vec- 
tor x{tk) of estimates of the drive system state, then 
the coder can calculate the best estimate x{tk+i) and 
choose Cfc-|_i = Cx{tk+i)- In this case the prediction error 
for a binary coder will be determined by the divergence 
rate of neighboring trajectories, i.e. relation H1U|) should 



be replaced by ^ exp{hTs) < A, where h > is the 
upper Lyapunov exponent of the chaotic drive system. 
This yields the bound < ln(2||C||)//i, instead of the 
bound Hill) . For the transmission rate it gives the nec- 
essary condition R* > h/ \n{2\\C\\) instead of the lower 
bound R* > Ly/A following from Hll() . If the condition 
R* > h/ \n{2\\C\\) holds, then the upper bound for trans- 
mission error A will decrease at each sampling interval 
[tk,tk+i) in ln(2||C||)) times and, therefore, will 

converge to zero exponentially. 

V. EVALUATION OF SYNCHRONIZATION 
ERROR 

Now let us evaluate the total guaranteed synchroniza- 
tion error Q — sup lim ||e(i)|| where sup is taken over the 

t — >oo 

set of transmission errors 5y (t) not exceeding the level A 
in absolute value The ratio Ce = Q/A (the relative er- 
ror) can be interpreted as the norm of the transformation 
from the input function 6y{-) to the output function e(-) 
generated by the system We will assume that the 
nonlinearity is Lipschitz continuous along all the trajec- 
tories starting from fi. More precisely, we assume that 

My)-ipiy')\\ <L^\y-y'\ 

for y = Cx, y' = Cx', dist(a;, fl) < A, dist(a;', ft) < A. 
The error equation (O can be represented as 

e^ALe + Cit), (20) 

where ||<^(i)|| < {L^ + ||L||)A. Choose L such that Al is 
a Hurwitz (stable) matrix and choose a positive-definite 
matrix P = > satisfying the modified Lyapunov 
equation PAl + AJ^P < —^P, for some > 0. After 
simple algebra we obtain the differential inequality for 
the function V{t) ^ e{tyPe{t): 

V < -iiV + e^P^(i) < -iiV + W ■ y/C^- 

Since V <0 within the set Vv>^i^^ sup y/Wf^PW)^ 

t 

the value of lim sup V{t) cannot exceed A^ (L^p + 

\\L\\) Aniax(-P)/M^- III view of positivity of P, 
\mUP)\WW < V{t), where Xmin{P), Xn^^P) are 
minimum and maximum eigenvalues of P, respectively. 
Hence 

Th^r||e(i)|| <C+A, (21) 

-here = vSff^- 

The relation 1)21(1 shows that the total synchronization 
error is proportional to the upper bound of transmission 
error A, i.e. can be made arbitrarily small for sufficiently 
large transmission rate R. 

One can pose the following problem: choose an optimal 
gain vector L providing the minimum value of Cg. How- 
ever an analytical solution is difficult to obtain in view 
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of the system nonlinearity. An alternative approach is to 
evaluate upper and lower bounds for Ce based on worst 
case inputs 6y (t) . Such a problem is similar to the energy 
control problem for systems with dissipation "2^, ^T?\ and 
Ce can be interpreted as excitability index of the system. 
Employing the lower bound for excitability index for pas- 
sive systems we conclude that if the gain vector 
L is chosen to ensure strict passivity of the system Q 
then the lower bound for Ce is positive, i.e. 



sup lim ||e(i)|| > C'A. 

|5„(t)|<A*^°° 



(22) 



Therefore for finite channel capacity the guaranteed syn- 
chronization error does not reduce to zero being of the 
same order of magnitude as the transmission error. Let 
us apply the above results to synchronization of two 
chaotic Chua systems coupled via a channel with limited 
capacity. 



VI. SYNCHRONIZATION OF CHAOTIC CHUA 
SYSTEMS 

System Equations. Consider the chaotic Chua system 
model: 



ii = p{-xi+ ip{y) + X2), t>0, 
±2 = xi - X2 + X3 (23) 
^X3 = -qx2, 
y{t) = xi{t), 



where y(t) is the sensor output (to be transmitted over 
the communication channel) , p, q are known plant model 
parameters, x = [xi, 0:2, ^3]^ g R'^ is the plant state vec- 
tor, the initial condition vector xo = x(0) is assumed to 
be unknown, ip{y) is a piecewise-linear function, having 
the following form: 

ip{y) = moy +\x + l\-\x-l\ 

+ 0.5(mi-mo)(|a; + l|-|x-l|), (24) 

where ttiq, mi are given plant parameters. 

Observer design. To obtain estimates x{t) of the cur- 
rent state x{t) of the system H23|l . the special case of a 
continuous-time observer jSJl is designed as follows 



(25) 



xi = p{-xi + (p{y) + X2) + heit), 
X2 = xi- X2 + £3 + heit), 
£3 = -qx2 + heit), 

[eit) ^ m - m, 

y{t)=xi{t), x(0) = Xo, 



where h, h, h ^re observer parameters, forming the 3x1 
observer matrix gain L = [h,h-,h]'^ ■ 



Subtracting from (g^) yields 

' ei ^ pi-ei + 62) + kiSyit) ~ ei) + ^lit), 

62 = ei - 62 + 63 4- h{5yit) " 62), 

63 = -962 + h {Sy it) - 63) , 

^Ut)=v{yit)-Syit))-^{yit)). 



(26) 



Equation 126|l describes the linear time-invariant (LTI) 
system e(t) = Ae(i), e(0) = xo — xo with the following 
matrix A: 



A = 



-p-h P 

1-^2 -11 
-/3 q 



(27) 



Matrix L should be chosen so that the observer H25|) 
stability conditions are satisfied, i.e. the characteristic 
polynomial -Dl(s) = det(sl — A^) is Hurwitz. For the 
observer (|25|l . the polynomial -Dl(s) has the form: 

£)l(s) = s^ + (1 + p + /i)s^ + i-q+ph + h)s 

-pq + hp-hq. (28) 

Evidently, we may find the matrix L for any arbitrarily 
assigned parameters di , d2, d^ so that the characteristic 
polynomial -Dl(s) — s'^ + dis"^ + d2S + d^. This leads to 
asymptotic convergence of the synchronization error e(t) 
to zero with prescribed dynamics in the disturbance-free 
case. 

Simulation results. For simulation the following pa- 
rameter values of the Chua system model (|23|l were cho- 
sen: p = 10.0, q = 15.6, mo = 0.33, nii = 0.22. The 
system exhibits a chaotic behavior, see yit) in Fig.Q] 

In our simulations parameter A has been taken from 
the set A = {0.1,0.2,0.5,0.7,1}. The sampling time Ts 
for each A has been chosen in accordance with (|17|l for 
Ly = 30 s^^. In the initial value Mq has been taken 
as Mq = 5, decay parameter p = exp(— O.lTs), and limit 
value Moo = M* = A/2. 

To evaluate the minimal synchronization error the op- 
timal observer gain matrix L*(A) was found numer- 
ically for several values of the transmission error A. 
We obtained L*(0.1) = [-4.66,0.50,-4.40]^, L*(0.5) = 
Ht-40, 0.46,-4.54]^ L*(1.0) = ^.97,0.46,^.47]^. For 
comparison the observer design by assigning a Butter- 
worth distribution of the observer matrix A^ eigenval- 
ues was performed. For the third order system the But- 
terworth design provides characteristic polynomial H28() 
as D(s) = s'^ -I- 2wos^ + 2uJqS + ujq, where parameter 
ciJo > specifies the desired estimation rate. In our ex- 
ample 070 = 6 s^^ is taken. It provides the observer 
eigenvalues: si = —6.0, S2,3 = —3.0 ±5.2 i. The observer 
feedback gain matrix L is found as L = [1.00, 5.54, 4.44]^^. 
For simulation the initial condition vectors for the sys- 
tems and were taken as xq = [0.3, 0.3, 0.3]'^ and 
io = [0,0,0]^. 

Simulation results for the coder (01, {Tj), lO and ob- 
server with optimally chosen gains for A = 1 are shown 
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FIG. 1: Outputs y{t), y{t), y{tk) time histories; A = 1. 




FIG. 2: Zooming of Fig[T]for t G [20, 22] s. 





FIG. 3: Time histories for A = 1: a) X2[t) ( ), X2{t) (— ); 

6) X,{t) (---), X3{t) (-). 



in Figs, n El The sampling interval is Ts = 0.02 s, 
which corresponds to the transmission rate J? = 50 bits 
per second. The following coder parameters were chosen: 
Mo = 5.0, Moo = 0.5, p = 0.998. It is seen that the syn- 
chronization process possesses sufficiently fast dynamics 
even in the presence of information constraints. 

It is seen from Fig. 0] that the difference in the limit 
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FIG. 4: Synchronization error Q vs A for different L. 




FIG. 5: Relative synchronization error Q/A vs A for different 
L. 
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FIG. 6: Synchronization error Q vs transmission rate 7?. 



synchronization error for different rationally chosen ob- 
server gains is not significant. Moreover, it is seen from 
Fig. [S] that the relative error does not approach zero for 
all choices of the observer gains. 

Dependence of the synchronization error Q on the 
transmission rate R is shown in Fig. demonstrating 
that the synchronization error becomes small for suffi- 
ciently large transmission rates. 
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VII. CONCLUSIONS 

We have studied dependence of the synchronization er- 
ror in the observer-based synchronization system both 
analytically and numerically. It is shown that upper and 
lower bounds for limit synchronization error depend lin- 
early on the transmission error which, in turn, is pro- 
portional to the driving signal rate and inversely propor- 
tional to the transmission rate. Though these results are 
obtained for a special type of coder, it reflects peculiar- 
ity of the synchronization problem as a nonequilibrium 
dynamical problem. On the contrary, the stabilization 
problem considered previously in the literature on con- 



trol under information constraints belongs to a class of 
equilibrium problems. As an intermediate result we ob- 
tained relation H19(l playing the role of an uncertainty re- 
lation between the transmission rate of information and 
the transmission error. 

Future research is aimed at analysis of controlled syn- 
chronization and control of chaos problems under infor- 
mation constraints. 

The work was supported by NICTA, University of 
Melbourne and Russian Foundation for Basic Research 
(project RFBR 05-01-00869). 
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